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Dynamical criticality has been shown to enhance information processing in dynamical systems, 
and there is evidence for self-organized criticality in neural networks. A plausible mechanism for such 
self-organization is activity dependent synaptic plasticity. Here, we model neurons as discrete-state 
nodes on an adaptive network following stochastic dynamics. At a threshold connectivity, this system 
undergoes a dynamical phase transition at which persistent activity sets in. In a low dimensional 
representation of the macroscopic dynamics, this corresponds to a transcritical bifurcation. We show 
analytically that adding activity dependent rewiring rules, inspired by homeostatic plasticity, leads 
to the emergence of an attractive steady state at criticality and present numerical evidence for the 
system's evolution to such a state. 
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Information processing systems often exhibit optimal 
computational capabilities when their parameters are 
tuned to critical states associated with phase transitions 
[iHl]. It therefore appears likely that our brains operate 
at criticality [5]. Although still hotly debated in neuro- 
science, the hypothesis of neural criticality is supported 
by recent experiments. Power-law distributions indica- 
tive of critical behaviour were observed in slices of rat 
cortex [SHE] as weU as EEG [SHTT], fMRI [Tg, and EcoG 
[llj measurements in humans. 

In the light of the experimental corroboration of neural 
criticality it is interesting to ask how a biological system 
can robustly self-tune its parameters to a critical state. A 
likely answer is found in the study of adaptive networks 
[131 [T3] , a class of models in which the dynamics on a 
network coevolves with the network structure. Already 
in 1998, it was noted that adaptive networks with slowly 
evolving topology can self-organize to a state where the 
dynamics on the network are critical [TS] . Adaptive self- 
organized criticality (aSOC) was subsequently demon- 
strated conclusively in a simple Boolean network model 
[IB] and then studied in detail in neural models [T7H23]. 

We note that insights in aSOC may not only advance 
our understanding of natural neural networks, but may 
reveal an important design principle for electronic com- 
puters: Within the next ten years, current manufacturing 
processes will hit fundamental boundaries [24j . Contin- 
ued progress will require utilising nanoscale components 
(e.g. nanotubes, nanowires, biomolecules) that cannot 
be positioned precisely with present-day techniques for 
manufacturing large-scale integrated systems. It is thus 
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conceivable that future computers may consist of active 
(nano-)clcmcnts that are deposited randomly and then 
left to self-tune to a critical state, where meaningful in- 
formation processing is possible. A promising mechanism 
for such self-tuning is provided by aSOC. Importantly, 
this mechanism does not require the rewiring of physical 
interconnections, but can be achieved already by local 
changes of conductivity between elements [22] that have 
recently been demonstrated experimentally |25| . 



Despite the prominent role of aSOC for information 
processing systems in both biology and technology, our 
understanding of the phenomenon is still limited. In 
many aSOC models, the critical state is identified by 
showing that certain quantities follow power-law distri- 
butions. However, power laws can appear due to other 
mechanisms besides criticality. It is thus desirable to 
make the corresponding dynamical phase transition di- 
rectly accessible to analytical investigation. 



The aim of the present paper is to develop a better 
understanding of aSOC by means of a simple conceptual 
model. We start in Sec. U with a brief review of the bi- 
ological background. In Sec. Illj we introduce a simple 
neuron model. In Sec. jIIH we use a moment closure ap- 
proximation ^26, for deriving a low-dimensional system of 
ordinary differential equations (ODEs) that capture the 
emergent dynamics. By analysing the bifurcation struc- 
ture of these ODEs, we show in Sec. |IV| that the model 
exhibits a non-equilibrium phase-transition. In Sec. |VJ 
we discuss the conditions under which the chosen topo- 
logical update rules drive the system toward the phase 
transition. Finally, in Sec. |VI| we show analytically and 
numerically that the model indeed exhibits aSOC. 



I. BIOLOGICAL BACKGROUND 

In biology, transmission of information between neu- 
rons occurs via cell contacts known as synapses. Up to 
the order of 10^ such connections (with different part- 
ners) can exist on a single neuron. Synapses allow a 
pre-synaptic neuron to depolarize a post-synaptic neu- 
ron, similar to polar devices that transmit current only in 
one direction. The topology of interconnections can thus 
be captured by a directed network, in which the nodes 
correspond to neurons and the directed links correspond 
to synapses. 

On a short timescale, neurons encode information in an 
electric potential across their cell membrane. In the ab- 
sence of inputs a neuron approaches a resting state with 
a characteristic membrane voltage. Depolarization of the 
membrane due to input from other neurons can lead to a 
firing state in which active mechanisms are used to emit 
a strong voltage pulse, which in turn excites connected 
neurons. After firing, the neuron enters a refractory state 
in which no excitation is possible before it finally returns 
to the resting state. 

On a longer timescale, the strength of synapses changes 
depending on the activity of the connected neurons. 
Thus, from an engineering perspective the synapse is 
a memristive element |27j . In biological terms the pro- 
cesses affecting synaptic strength are collectively known 
as synaptic plasticity. Here, we specifically consider 
homeostatic synaptic plasticity (HSP) [281129] . which de- 
creases the strength of synapses if the activity of a neuron 
is high, and increases the strength if the activity is low. 



II. DISCRETE NEURAL MODEL 

In the present paper we consider a directed network of 
N nodes. At any time a given node is in one of three 
different states: resting (inactive, I), firing (F), or refrac- 
tory (R). We start from a random network in which a 
majority of the nodes are in the inactive state, whereas 
a small fraction is in the firing state. The node states 
are then evolved according to the following rules: Firing 
nodes become refractory at a rate i, refractory nodes be- 
come inactive at rate r, and, for every link pointing from 
a firing neuron to an inactive neuron, the target neuron 
is set to the firing state at rate p. 

The network topology is evolved according to a rule 
modelling synaptic plasticity: A firing node looses an 
incoming link at the rate / whereas new links are estab- 
lished between nodes at the rate g. The creation and 
deletion of links is reminiscent of the formation of synap- 
tic contacts in the developing brain, but is used here as a 
discretized model of the continuous changes of synaptic 
weight in the adult brain. By formulating the model in 
terms of discrete linking and unlinking events we avoid 
additional complications caused by real-valued link dy- 
namics. 

The conceptual model described above is one of the 



simplest conceivable settings that allows studying the in- 
terplay between the spreading of excitation and home- 
ostatic topological evolution. Excitation dynamics are 
modelled by stochastic transitions with constant rates. 
As a consequence, a node can be excited by a single ac- 
tive neighbour, whereas in real neural systems, multiple 
inputs are needed. Moreover, the time a node spends in 
firing state is exponentially distributed, while real action 
potentials have a well defined length. However, as known 
from epidemiological modelling, none of these simplifica- 
tions impair the validity of the model predictions [30,. 



III. LOW-DIMENSIONAL APPROXIMATION 

Let us now derive a low-dimensional description for the 
time evolution of macroscopic quantities in the limit of 
infinite N. To this end, we employ a moment closure ap- 
proximation (MCA) [3^1 EH [32j ■ We denote the densities 
of nodes in a certain state by [F], [I], [R], respectively. 
Similarly, we denote the per-neuron density of links from 
a node in state X to a node in state Y by [XY] . 

Consider [F], the density of nodes in firing state: It 
increases when a firing neuron excites an inactive node, 
which happens at rate at rate p[FI], and decreases when 
a firing node becomes refractory, which happens at rate 
i[F]. By such reasoning, we obtain equations for the time 
evolution of [F], [I] and [R], which however, do not form 
a closed system as they depend on the link density [FI]. 

One possibility to close the system of equations is to 
approximate the link density by the densities of the nodes 
involved, usually by assuming [FI] « fc[_F][/], where k 
is the mean degree of the network. Note that in such 
a m,ean-field approximation, correlations between neigh- 
bouring nodes are neglected. 

Alternatively, it is possible to treat the link densities 
themselves as dynamical variables and derive evolution 
equations for them. In general, these equations depend 
on triplet densities, which can then be approximated by 
link and node densities. This approach is often referred 
to as the pair approximation; it is described in detail in 
Appendix IX] 

Using the pair approximation and assuming a Poisso- 
nian degree distribution, we obtain the following ODE 
description of the system: 



[F] = p[FI] - i[F] 
[R] = i[F] - r[R] 
[FF] = -2^[FF]+p(iI^ 
+ .9[^1[^1 
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r[RF] - l[RF] 



[RF] =i{[FF] 

+ 9[RW] 
[RI] = i[FI] -pimm+r{[RR] - [RI]) + g[R][I] 

[RR] = i {[FR] + [RF]) - 2r[RR] + g[R][R] 
k = g~l[F] 
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(Ij) 
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where the last equation exphcitly captures the change of 
the mean degree k of the network. Writing the equation 
for k saves us from writing the longer equation for IF, 
due to 



[IF] = k- [FF] - [FI] - [FR] - [II] - [IR] 
- [RF] - [RI] - [RR] . 



(2) 



Similarly [/] follows from the conservation relation for 
nodes 



[I] = l-[F]-[R] 



(3) 



IV. NON-EQUILIBRIUM PHASE TRANSITION 

Before we address the dynamics of the full system, let 
us first focus on the case I = g = Q, where the network is 
static. In this case, the right hand side of the equation 
of motion for k vanishes such that k becomes a control 
parameter. Below, we will refer to the system without 
topological evolution as the static network model and to 
the complete system as the adaptive network model. 

The static network has a trivial steady state x'^ = 
([F]'',--- ,[RR\^) in which all nodes are inactive, i.e 
[If = 1, [//]0 = fc, [F]0 = [i?]0 = [i^F]0 = . . . = [i?i?]0 = 
0. Depending on parameters, it may moreover have a 
non-trivial, active steady state, in which [/] < 1. The 
stability of the trivial steady state is determined by the 



spectrum of the Jacobian matrix J|xO G 



plOxlO 



, where 



Jij = dxi/dxj. The steady state is asymptotically stable 
if all eigenvalues of J|xO have a negative realpart 33!. 

If the variables Xi are ordered as in Eqns. (fl]), the non- 
vanishing entries of Jj^o are 
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Jl,l = J7,3 ~ JtA = "^7,5 


= Jrfi ~ 


--J7,s 


JlA = "^3,4 = P , 






^2,1 = ^5,3 = Js,3 = J9,4 


= ^^10,5 


~ Jio 


•^2,2 = "^4,5 ~ J(>,7 ~ Jg,9 


= •^9,10 


= r , 


J'i/i = -2i , 






^4,4 = -i + (fc - l)p , 






•^5,5 = Jt,? ~ Js,8 — —i - 


r , 




J6,4 = -2kp , 






Jt.io = —i + r , 






J9,9 = -r , 






iio.io = — 2r . 







In the following, we assume i > 0, r > and p > 0. 
The characteristic polynomial can then be factored into 



7 linear factors with negative real roots and a remaining 
third order polynomial 

P{X) = (-r-A)3(-z-A)2(-i-r-A)(-2r-A)/(A) , (4) 

where 

/(A) = I (t^ + 2i ((1 - k)p + r) + (1 - 2fc)pr) 

+ {hi^ + (1 - fc)pr + ii ((1 - kp) + r)) A (5) 
+ (4z + (1 - k)p + r) A^ + A^ . 

In order to assess the stability of x° without explicitly 
calculating the roots of /(A), we use the Routh-Hurwitz 
theorem [34 . It states that the roots of a polynomial 
p{x) = x" -I- 6ix"~^ + • • • + hn-\x + hn all have negative 
real parts if the Hurwitz determinants 



A. = 



bi 
b3 



1 
62 



>2k- 



^1 b 



2fc-2 





61 

b2k- 



'J2k-i 
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with k = 1 . . . n, are all positive. Calculating the Hurwitz 
determinants of the Jacobian matrix Jj^o and evaluating 
the positivity conditions reveals that the trivial steady 
state is stable if and only if 



k< 



r/2 



P 



(7) 



At fc = fcc the trivial steady state becomes unstable as 
the systems undergoes a transcritical bifurcation. In this 
bifurcation a nontrivial steady state enters the positive 
cone of the state space, becoming a physical solution. 
The transcritical bifurcation thus marks a transition be- 
tween the trivial inactive state and an active state in 
which ongoing activity is observed. This transition is in 
the same universality class as directed percolation. 

The role of the transcritical bifurcation is illustrated 
in a representative bifurcation diagram in Fig. IT] The 
diagram shows that the analytically predicted bifurca- 
tion point is in very good agreement with numerical re- 
sults from agent-based simulations of the system. Fur- 
ther, numerical continuation of the ODEs (II]) shows 
good agreement close to the bifurcation point. By con- 
trast, closing the MCA at mean-field level yields [F]'^^ = 
r{l — i/{pk))/{i + r) for the non-trivial steady state. This 
can be seen to be a much poorer approximation and un- 
derestimates the bifurcation point. 



V. LOCAL INFORMATION AND TIME-SCALE 
SEPARATION 

The previous section showed that the static network 
model exhibits a transcritical bifurcation. In the lan- 
guage of statistical physics this bifurcation constitutes a 
phase transition. To establish that the adaptive network 



Et, 




FIG. 1. Bifurcation diagram for the static network. Plot- 
ted is the steady state density of firing neurons [F] over 
the network's mean degree k. The solid line marks stable 
steady states of the static system, the dashed line unstable 
ones. At k = kc ~ 5.6, the inactive steady state loses sta- 
bility in a transcritical bifurcation. The respective transi- 
tion from an inactive to an active phase is already observed 
in individual-based simulations with A'' = 10^ neurons (cir- 
cles). Note that the critical k is nicely predicted by the 
link- level approximation Eqns. ([l|, but underestimated by 
a MCA at mean- field level (dotted lines). Parameters used 
were p = 0.2, i = 0.95, r = 0.4. Nontrivial steady states were 
calculated using AUTO [35]. 



is argued that in adaptive networks robust self-organized 
criticality is possible because the dynamics on the net- 
works make information on the global topology available 
in every node. Based on this information the nodes can 
then infer the global phase and adjust the local topology 
accordingly. Indeed, in previous publications [TBI HHI the 
networks nodes extracted the global phase information by 
integrating their dynamics over a long time. While some 
biological processes are known to enable such temporal 
averaging |361 I37j , we explore an alternative explanation. 

In our model, the local accessibility of the global phase 
information is limited, as both the dynamics on the static 
network as well as the topological updates are memory- 
less processes that depend only on the current state of a 
single link or node. Hence, neurons in the resting state 
do not possess any information about the global phase, 
as they occur in both phases. By contrast, neurons in 
the firing state can infer the global phase from their lo- 
cal state, as their occurrence is restricted to the active 
phase. To achieve aSOC despite the limited local accessi- 
bility of the global phase, links have to be created tenta- 
tively as long as no definitive information is known, but 
destroyed decisively once information about the global 
phase is available. This corresponds to a separation of 
timescale between the link creation and link deletion pro- 
cess 



<1 . 



(9) 



model shows self-organized criticality we have to show 
that the evolution of the connectivity drives the system 
to the critical point. 

In the discussion leading up to Fig. [T] we treated k as 
a parameter of the system. For studying the evolution of 
the connectivity we now consider fc as a dynamical vari- 
able that evolves according to Eq. (Ik). By introducing 
dynamics in k we change the dynamical system, which 
can potentially lead to a changed bifurcation diagram. 
Indeed, in general a diagram analogous to Fig. [l] cannot 
be drawn for the adaptive network model because k is 
not available as a parameter axis anymore. 

The pitfall described above is inherent in the concept 
of SOC: Criticality can only be cleanly defined in a dif- 
ferent system, where the self-organization is absent. To 
circumvent this pitfall one has to demand that the self- 
organization acts on a slower timescale. In the example 
studied here this implies. 



g,l<^p,i,r 



(8) 



In this case, the bifurcation diagram of the static network 
model reappears as the bifurcation diagram of the fast- 
subsystem of the adaptive network model, where k can 
again be treated as a parameter. 

The genesis of aSOC in our model requires additionally 
a second time scale separation. To see this let us revisit 
a widely held plausibility argument for aS0C[T31[TS]. It 



Note, that the twofold separation of time scales consti- 
tuted by Eq. ^ and Eq. ^ is typical for SOC models 
[T8l [T9l [22] , although it is sometimes hidden in a non- 
local model definition [3S| |31] . Our analysis reveals that 
it is indeed a necessary ingredient to achieve criticality 
without centralized control. 



VI. SELF-ORGANIZED CRITICALITY 

Let us now show that the adaptive network model has 
a stable steady state, which in the double limit / -^ 0, 
e — > coincides with the bifurcation of the static network 
model. We start by defining x* = {[F]*,- ■ ■ , [RR]*,k*) 
as a steady state of the adaptive network model (IT]) . Eval- 
uating the stationarity conditions fc|x* — 0,[F]\x* = 
and [R]\x' = 0, we can immediately read off the steady 
state values 



[F]* = e, [R]* - et/r, [FI]* = ei/p 



(10) 



Therewith, the remaining equations ([i^F]|x* — 

0, • • • , [i?i?]|x* — 0) become linear and can be solved by 
a computer algebra system. Due to the time-scale sepa- 
ration, higher-order terms in e and I can be dropped, and 



we obtain 



[FFr = -^e 



[FRr = [RFr^^^e 
2 t + r 

iiRr = l(k. + l)e 



[Riy = -ikr. 

r 



[RR] 



1 i 



2 r i + r 
k* = kc + 



4i{i + r) 



r V 2 / I + r 



(11a) 
(lib) 
(lie) 
(lid) 
(lie) 
(llf) 

(llg) 



For e,l <^ 1, the steady state x* is always stable, which 
can be shown by calculating the characteristic polyno- 
mial of J|x* and then checking the Hurwitz determinants 
for positivity (again, dropping higher-order terms in e 
and I). Moreover, the Eqns. (11 1 reveal that in the double 
limit / ^ 0, e ^ 0, [//]* ^ fee, while [F]* , ..., [RR]* -^ 0, 
i.e., X* converges towards the transcritical bifurcation of 
the static network model. 

Recall that the ODE system describes the model 
system in the limit N —i' oo; as observed recently, this 
is also the limit in which criticality may be expected 
for non-conserving dynamics [lOl SI] • We thus conclude 
that for N — > oo, and ^ — >■ 0, e — ?► 0, the HSP-inspired up- 
date rule self-organizes the adaptive system to criticality. 

To confirm aSOC in the discussed model, we ran 
individual-based simulations. We start with a network 
of N nodes, each of which has an outgoing connection 
to any other node with probability ko/N. We initial- 
ize 5% of all nodes in the firing state, all others in the 
inactive state. Then, nodes and links are evolved accord- 
ing to the rules described in Sec. |Tl] using the Gillespie 
algorithm [12]. 

For finite N, the simulations tend toward the absorbing 
inactive state due to demographic stochasticity. To com- 
pensate for the finite size effect, we include an additional 
process: Inactive nodes fire spontaneously at rate s. This 
process has an immediate biological interpretation as it 
reminds of the spontaneous activity observed in neurons. 
To reconcile the simulations with the low-dimensional de- 
scription (fTl), s has to vanish in the N ^ oo limit. A 
plausible assumption is s = c/N, where c can be chosen 
arbitrarily. We have verified that for sufficiently large 
system sizes the particular choice of c does not signifi- 
cantly influence the evolving connectivity (cf. Fig. [2]). 

We start in Fig. [3] by plotting the time evolution of the 
mean degree k in networks with different initial config- 
urations. All networks approach the same connectivity, 
which corroborates that the adaptive network model has 
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FIG. 2. (Colour online) Mean degree k of evolved networks 
for different network sizes A'' and different spontaneous firing 
rates s = l/(1000Ai') (yellow triangles pointing right), s = 
l/(100Af) (black triangles pointing left), s = 1/{10N) (blue 
squares), s = 1/7V (green circles), s — 10/ N (red triangles 
pointing up). The dotted line marks the analytical steady 
state value for I — 0.01, e = 10~*. Simulations were run 
for lO'^ time units, plotted connectivities are averages over 



the last 5 ■ lO"" 
0.95, r = 0.4. 



time units. Other parameters: p — 0.7, i 



■i^ 3 




FIG. 3. Time evolution of the mean degree of networks 
(TV = 10000,/ = 10"^e = 0.01, s = 10"*), starting from 
different initial connectivities. Inset: In-degree distribution 
of an evolved network after 5 • 10® time steps {N = 10^,1 = 
10~'',e — 10~^,s — 10~^, open circles) and Poissonian distri- 
bution around the same mean (filled circles). Other parame- 
ters: p = 0.7, i = 0.95, r = 0.4 in both cases. 



exactly one stable steady state to which it converges ir- 
respective of initial conditions. 

The inset of Fig. [3] shows that the in-degree distri- 
bution stays Poissonian throughout topological evolu- 
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agrees with the analytical estimate of the critical state 
except for a small but systematic deviation. This dis- 
crepancy can be understood considering Eq. (llgl. The 



FIG. 4. Mean degree k of evolved networks for different values 
of p. The dashed line marks the analytically obtained criti- 
cal connectivity kc- Simulations were run for 10* time units. 
Parameters: iV = 10^,1 = 10"®, e = 0.0015, s = lO"'',* = 
0.95, r = 0.4 
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FIG. 5. Average density of firing neurons (after transients) 
around the evolved connectivity (marked by vertical dashed 
line). Values above or below the evolved connectivity were 
obtained by adding (removing) random links to (from) the 
evolved topology. The network was evolved over 5 ■ 10® time 
units. Parameters: N = 10^, p = 0.7, i = 0.95, r = 0.4,/ = 
0.01,6 = 10"*, s = l/(100iV). 



tion and thus retrospectively corroborates the assump- 
tion made in the derivation of the MCA. 

So far we have shown that the network self-organizes 
to a unique value of the connectivity. It remains to test 
whether this value is indeed the critical connectivity. As 
a first test, we compare the connectivity of the networks 
evolved in numerical simulations with the critical con- 
nectivity kc predicted by the analytical approximation 
from Eq. ([7]) . As shown in Fig. m the numerical result 



second and third term of this equation are positive for all 
i,r,kc,l,e > 0. Hence, the systematic deviation k* > kc 
is consequence of the networks' topologies being evolved 
with small but finite rates e, /. 

As a second test, we directly probe the dynamics on 
the evolved networks. For this purpose, we first evolve 
the topology until the mean degree has reached a sta- 
tionary state. Then, we deactivate the adaptive addition 
and deletion of links but manually add (delete) links from 
the network. After every addition (deletion) we let the 
dynamics on the network reach a stationary level and 
record the average number of active neurons. As shown 
in Fig. [5] this procedure recreates the phase diagram of 
the system numerically. It thereby provides direct evi- 
dence of the criticality of the evolved state. The slight 
displacement from the critical point can be understood 
recalling the analytical results given above: According 
to Eq. (10), [F]* = e in the evolved network. Hence, the 



displacement toward the active regime can be attributed 
to the small but finite rates e, I used in the simulations. 



VII. DISCUSSION 

In summary, we have shown that activity-dependent 
synaptic plasticity self-organizes a neural network to crit- 
icality, provided that driving is slow and there is an ap- 
propriate separation of timescales between potentiation 
and depression. 

We have considered a simple discrete neural network 
model and used moment-closure approximation to derive 
a low-dimensional ODE representation, in which the dy- 
namical phase transition becomes manifest as a transcrit- 
ical bifurcation. By adding activity dependent rewiring 
rules we obtained an adaptive network model that has 
one attractive steady state. We have shown that, in the 
limit of infinitesimally slow topological adaptation, this 
steady state coincides with the bifurcation of the static 
network model, and, thus, that the adaptive model dis- 
plays aSOC. 

We emphasize that the particular choice of the param- 
eters p, r, and i does not affect the qualitative model 
behaviour but only the quantitative predictions, in par- 
ticular the evolved connectivity: The probability that a 
firing node excites a specific neighbouring node is p/i. It 
is apparent from Eq. [7] that choosing higher values for 
p/i will result in lower evolved connectivities kc- In or- 
der to keep simulations computationally feasible, which 
specifically means maintaining a low number of links, we 
have chosen unphysiologically high values for p/i. It can 
be seen from our analytical treatment that this affects 
neither the existence of a steady state at criticality nor 
the stability of this state. 

Our work provides a conceptual angle on self-organized 
criticality of neurally inspired models. Although the phe- 



nomenon has been demonstrated in several earlier works 
[m |43] , the simplified model studied here is the first sys- 
tem in which it is analytically tractable. Using dynamical 
systems and bifurcation theory, aSOC can be established 
more rigorously than in more realistic models and exper- 
iments, where evidence for criticality is mainly provided 
by the numerical observation of power laws. We believe 
that this approach will prove useful in the context of 
other, more complex dynamical phase transitions. For 
example, it is conceivable that a similar strategy may 
help to elucidate the recently observed self-organization 
to the onset of synchronous activity [21] ■ We hope that 
the approach can thus contribute to settling an ongoing 
debate regarding the existence and function of critical- 
ity in the brain [5] and potentially to the development of 
self-organising electronic circuits. 



Appendix A: Moment-closure approximation 



As described in Sec. Ill the time evolution of node 
densities is given by 



[F]=p[FI]-t[F], 
[R] =i[F] -r[R]. 



(Ala) 
(Alb) 



For an expansion beyond mean-field level, we need to 
derive equations for the time evolution of link densities. 
There are four ways a link of a certain type be created 
or destroyed: 1) One of the link's nodes changes inde- 
pendently of others, 2) the link is added or removed as 
a result of topological evolution 3) one of the nodes is 
excited by the other via the link in question, or 4) one of 
the nodes is excited via a different link. Processes of type 
1), 2) and 3) can be understood on the level of nodes and 
links. To understand 4), consider that activity in a given 
focal link, say, excitation along an FI link, does not only 
affect the links itself, but also links connecting to it, here 
for instance other links connected to the I-node. Thus, 
processes of type 4) can only be captured if the density 
of subgraphs with two links, so-called triplets are taken 
into account. 

Consider the change in the density of FI links. 



[FI] = p{[F>I>I] - [FI] - [F>I<F] 
^i[FI]+r[FR]+g[F][I], 



(A2) 



where [X>Y>Z] denotes the density of directed triplets, 
with the directionality of links indicated by "<"and ">" , 



respectively. The three rightmost terms are of type 1) 
and 2). The others can be understood as follows: a 
F>I >I triplet turns into a F>F >I triplet at rate p, in- 
creasing the number of FI links by one (the underline 
denotes the link on which the update rule is applied). 
At the same rate, the F in FI links excites the I, which 
decreases the number of FI links by one. Finally, F>I <F 
triplets change to F>F <F at rate p, destroying the FI 
link not underlined. In the last transition, there are ac- 
tually two FI links being destroyed, but we have already 
counted one of them as a process of type 3). Evolution 
equations for other types of links are obtained using the 
same reasoning. 

The pair approximation closes the system of equations 
by approximating the occurring triplets densities in terns 
of link densities. Let us start with a triplet of the type 
X>Y>Z. It consists of one XY link, which we know oc- 
curs with densities [XY] . If we assume that the states of 
next-nearest neighbours are uncorrelated, we can calcu- 
late the expected number of links from the Y node to a Z 
node as [yZ]/[F]. We can thus approximate the triplet 
density as 



[X>Y>Z] 



[XY] [YZ] 
[Y] ■ 



(A3) 



Triplets of the type X>Y<Z are approximated simi- 
larly: Again, the triplet contains an XY link, which oc- 
curs with density [XY]. To calculate the expected num- 
ber of ZY links connected to the central Y node, we now 
need to take into account that we already know that one 
incoming link is a XY. Thus, we need to consider the 
mean excess degree q, i.e. the expected number of in- 
coming links that the Y has in addition to the XY link. 
Assuming that each of these links is a ZY link with prob- 
ability [Zy]/(fc[F]), we can write 



[X>Y<Z] 



q[XYWZ] 
k [Y] ■ 



(A4) 



Note that for symmetric triplets (X>Y<X), Eq. (A4) 
has to be modified to avoid double-counting: 



[X>Y<X] 



q_[XYl 
2k [Y] 



(A5) 



Assuming that the in-degree distribution of the evolv- 
ing network is Poissonian, we can simplify the triplet 
approximations further. For Poissonian degree distribu- 
tions, q = k [53|, which leads to the expressions in Eq. (II]). 
Note that the assumption is confirmed by the numerical 
results in Sec. 
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